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Colloidal gel aging is investigated using very long duration brownian dynamics simulations. The 
Asakura Oosawa description of the depletion interaction is used to model a simple colloid polymer 
mixture. Several regimes are identified during gel formation. The Intermediate scattering function 
displays a double decay characteristic of systems where some kinetic processes are frozen. The 
(3 relaxation at short times is explained in terms of the Krall-Weitz model for the decorelation 
due to the elastic modes present. The a relaxation at long times is well described by a stretched 
exponential, showing a wide spectrum of relaxation times for which the q dependence is r a = q~ 2 ' 2 , 
lower than for diffusion. For the shortest waiting times, a combination of two stretched exponentials 
is used, suggesting a bimodal distribution. The extracted relaxation times vary with waiting time as 
r a = t°' 66 , slower than the simple aging case. The real space displacements are found to be strongly 
non-Gaussian, correlated in space and time. We were unable to find clear evidence that the gel aging 
was driven by internal stresses. Rather, we hypothesise that in this case of weakly interacting gels, 
the aging behaviour arises due to the thermal diffusion of strands, constrained by the percolating 
network which ruptures discontinuously. Although the mechanisms differ, the similarity of some of 
the results with the aging of glasses is striking. 

PACS numbers: 82.70.Dd,61.43.Hv,61.20.Lc 



I. INTRODUCTION 

Like other non-equilibrium phenomena, colloidal gela- 
tion has attracted considerable interest in the last twenty 
years. Colloids are found in many forms and are highly 
relevant for industry. Muds, paints, cosmetics, and many 
food products are examples where understanding of col- 
loid science is essential in the design of product proper- 
ties. Studies can be conducted over a wide range of den- 
sities, and the interactions can be finely tuned in terms 
of range and depth of potential by altering the particle 
coating, solvent conditions, or solutes present. 

Importantly, under certain conditions it is possible to 
make colloid particles form a macroscopic gel. Although 
there is no universal definition, a gel is usually under- 
stood to be a percolating viscoelastic network that has 
no long range order (although it may have some local or- 
der), but a characteristic length scale is present; gels are 
usually fractals at low volume fractions but this is not 
a necessity. Common examples of colloid gels include 
yogurt, toothpaste, and clays. 

The system under study here is a colloid polymer mix- 
ture. Such systems have been known to have interest- 
ing equilibrium and non-equilibrium properties for many 
years [Q ; they are also a stepping stone to understanding 
the more complicated mixtures that make up many ordi- 
nary products: shampoo and paint are examples which 
involve complex polymer colloid mixtures. Such systems 
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can display depletion flocculation, and the resulting gels 
are often used as one of the simplest models of weakly 
attractive system due to their relative simplicity arising 
from the absence of charge effects. 

The properties of non-equilibrium systems as they age 
are not only relevant for applications, but are also a cru- 
cial step to understanding the very nature of gels. Indeed 
if slow aging were not a feature of the gels, phase sepa- 
rations would only result in equilibrium phases. 

The dramatic slow down of the dynamics in gels has 
prompted investigations into their possible link with 
glasses. While gels are usually found at lower volume 
fractions and when strong interparticle interactions are 
present, colloid glasses are usually found at high volume 
fractions, above <f>— 0.58, or for deep quenches. They 
share many properties such as being out of equilibrium, 
non-ergodic, the absence of long range order and slow 
dynamics. Since work on the glasses and their aging is 
extremely topical in soft condensed matter, this work will 
look for similarities between the aging of colloid gels and 
glasses. 

There is abundant literature on gels and on colloid 
polymer systems in particular. Asakura and Oosawa 
were the first to describe the attraction between colloid 
particles when non-adsorbing polymer is present in solu- 
tion The depletion mechanism which is responsible 
for this can be understood in different ways. The essen- 
tial point is the presence of a region around the colloids 
from where the polymer coil centre of mass is excluded. 
When two colloid particles come together the regions of 
excluded volume overlap resulting in an increase of the 
free volume available to the polymer, causing an increase 
of their entropy and an overall decrease of the free en- 
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ergy of the system. The problem can also be analysed in 
terms of the net osmotic pressure on a pair of colloids as 
they come together. 

Theoretical work on the system was done by Vrij ||, 
who independently rediscovered the Asakura and Oosawa 
model, and demonstrated the presence of spinodal insta- 
bility in the system. Gast et al. treated the polymer as 
a perturbation to the hard sphere results thus producing 
the first phase diagrams j|, |[. Lekkerkerker et al. || 
used a free volume approach to treat this problem while 
Djikstra et al. formally integrated out polymer degrees 
of freedom 0. It has been shown that for size ratios 
below k = 0.154 there can be no three body interactions 
and that the two body Asakura and Oosawa (AO) poten- 
tial can therefore be an adequate description. A review 
of the details of the different approaches to the colloid 
polymer system can be found in [Q. It is important to 
note that the depletion potential give well controlled sys- 
tems where the interparticle interaction potential well 
is usually in the range of 1 — 10/cT, very different from 
cmultions, for example, where interaction energy can be 
of order lOOfcT. 

Felicity et al. || and Lodge et al. || modelled col- 
loid gels interacting via a Leonard- Jon es p otential and 
studied their viscoelastic properties || |HJ; Bijsterbosh 
et al. investigated fractal colloidal gels for various in- 
teractions using brownian dynamics simulation methods 
pl| . Dickenson et al. |L2||, and Bos et al. |l3j have 
worked on stress relaxations after shear stresses were ap- 
plied. Soga et al. have worked on percolation jl4|] and 
phase transitions |l5| in the colloid polymer system. The 
yielding behaviour of 2D colloid gels was investigated by 
West et al. (l6 . 



Experimentally, most notable was the work of Grant 
and Russcl |r^| , Poon et al. M, Lekkerkerker et al. [ 18[ , 
and Verdin and Dhont]l9|]. An important result from 
all these works was that at moderate volume fractions 
(0.1 < (f> < 0.4) gelation can be expected to occur for 
high polymer concentrations (deep quenches). Also, for 
ranges of interaction of less than a third of the particle 
diameter no liquid phase formed; the liquid gas binodal 
becoming metastable with respect to fluid solid phase 
separation. In light scattering, the gelation is marked 
by the 'freezing' of the low q peaks demonstrating that 
the dynamics become arrested as the sol gel transition is 
reached. 

The following experiments are explained in more detail 
since they contain the most interesting data for compar- 
ison with our simulations. Cipelletti et al. pfj|] have pro- 
duced experimental results on the aging of polystyrene 
colloid gels at low volume fractions (10~ 4 < cj) < 10 -3 ). 
The gels where formed by aggregation due to the van 
der Walls forces, and were found to be fractal. The ag- 
ing behaviour was sampled for waiting time going from 
hours to ten days, the relaxations where stretched with 
the exponent greater than one. The unusual q depen- 
dence of the relaxation times, Tf oc q , was explained in 
terms of a dipole elastic field caused by the high internal 



stresses of the gels. The waiting time dependence of the 
relaxation times was nearly linear, Tf oc t^ 9 , following 
a period of exponential growth (this exponential regime 
was also observed in an aging clay suspension pj| ). 

The same group, in a series of rheology and dynamic 
light scattering experiments, also found the same be- 
haviour in a system composed of multilamellar vesicles, 
where, once again, it was thought to be the strong inter- 
nal repulsive stresses that dominate the response |2l| . 

Experiments on the dynamics of colloidal gels have 
been performed by Krall and Weitz |^3| for the short 
time (or /3, see below) relations. They also presented a 
model for the elastic modes present in a fractal gel, from 
which the mean squared displacements of strands con- 
strained by the network can be calculated. Working on 

tat densies conparable with this work, Romer et al. 
, showed that the predictions of this model held even 
for non-fractal structures. 

Both short and long time scale dynamics where ob- 
served by Solomon et al. Q in their study of adhe- 
sive colloids at low volume fractions and for interparti- 
cle interactions strengths close to that used in this work 
( lOfcT). The Krall- Weitz model was used at short times, 
which a simple exponential was observed at long times. 

A commonly calculated quantity is the intermediate 
scattering function (ISD), a time dependent correlation 
function. When looking at glasses and gels the self part 
of the ISD will display a double decay. At short times a 
fast decay often referred to as the (3 relaxation, and at 
long times a second decay called the a relaxation. Struc- 
tural arrest, brought about by increased density for hard 
spheres or by quenching for a Leonard- Jones system, will 
lead to a wide separation of these time-scales and the ap- 
pearance of a plateau at intermediate times see reference 
p6| and references therein. Since the scattering function 
is wavevector dependent, spatial information can also be 
obtained. 

The a and (3 relaxations are often explained in terms 
of the cage effect. For a system of purely repulsive parti- 
cles, a particle is confined by its neighbours which form 
a cage around it; the particle is free to diffuse within its 
cage, thus setting the (3 time scale. The time required 
for the particle to diffuse out of the cage is the a time 
scale. As a system becomes denser, cage rearrangements 
become increasingly difficult, causing the a relaxations 
to occur on longer time scales. For systems of attrac- 
tive particles, one can talk about an open cage, whereby 
a particle would be confined to its neighbours because 
of the attractive forces; two relaxation time scales also 
result from this. On this basis it has been argued that 
the gel transition should be described in terms of Mode 
Coupling theory (MCT) as an ergodicity breaking tran- 



sition l27[ |28| |29[ |30|. MCT has considerably helped 
to understand supercooled liquids p6| , and high volume 
fraction colloid glasses j3lfl , and its predictions are now 
being tested for gels as well plfl . The kinetic arrest of 
gels has also been reported by Segre et al. , who com- 
pared it to the features of the hard sphere colloid glass 
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transition |3g |. 

Ideal MCT assumes closeness to equilibrium and is 
therefore only useful for the fluid-gel transition. Non- 
equilibrium MCT has also been developed by Bouchaud 
et al. p3| which predicts aging behaviour. Some of these 
predictions have been qualitatively tested on a Laponite 
glass pi]] . However, earlier theoritical models have inter- 
esting aging properties. Cugliandolo and Kurchan have 
solved the spin glass model for the out of equilibrum dy- 
namics J34|. Crucially, these equations are identical to 
those produced by some MCT models. 

Another approach to explain the behaviour of glasses 
has been the study of the properties of the Potential 
Energy Surface (PES) as proposed by Goldstein p5| . 
One model which uses this concept is the successful trap 
model of Bouchaud [fjrj which has been widely applied 
for spin glasses and of which many variants exist, this 
was one of first models to make useful predictions about 
aging. Another formalism is now introduced to facilitate 
discussion. The Inherent Structures (IS) of a system are 
defined as the local minima of the PES; all points in con- 
figuration space which end up in the same IS following a 
steepest decent path define the basin of that IS |37j. Sci- 
otino and Tartaglia distinguish between four dynamical 
regimes: a high temperature region where the system is 
fluid; an intermediate region where the system is always 
close to the ridges separating the different basins, this 
being the regime in which MCT is successful; the third 
regime is found at low T where the system populates IS 
with lower and lower energy. When the time of escape 
from a basin becomes the longer than the experimental 
time scale the system no longer equilibriates, which is the 
fourth regime, normally called a glass |38| . In this for- 
malism the two step relaxation of the correlation function 
is seen as decoupling of the fast (intra-well) vibrational 
excitations and the slower rearrangements corresponding 
to the system sampling different basins which for low T 
may require activated processes. 

It is useful to set up a few hypothesises about the possi- 
ble causes of the aging in gels. Firstly, very slow crystalli- 
sation about some locally highly ordered regions is a pos- 
sibility although this is not normally seen in experiments. 
Secondly, the aging may be stress driven as argued to be 
in some experimental systems |2^, p2| . Thirdly, thermal 
fluctuations and random breaking of bonds to allow the 
relaxation of the network as suggested by Solomon et al. 
p5j. Compactification of strands up to a point where 
they might dislocate has also been suggested as a possi- 
bility to explain the collapse of transient gels |39|. 



collisions of the solvent molecules with the particle. The 
polymer coils are not simulated, their presence only being 
felt through the attractive potential between the colloids. 
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Where M and X are the inertia and position matrices 
respectively; F c are the conservative forces, F B are the 
random brownian force, and F D are the dissipative drag 
force. This last term encompasses all the hydrodynamic 
interactions, which couples the forces on the particles 
through the solvent. It is these forces which oppose the 
relative motion of the particles, and are most important 
in flow. They will affect non-equilibrium phenomena at 
rest, but in common with many others we neglect the 
many body interactions in this study. In this model, this 
term is reduced to Stoke's drag for a spherical particle 
(a), applied individually: this is the free draining ap- 
proximation. So the drag force on each particle satisfies: 



fi D = -aVi = -Girrjavi 



(2) 



where rj is the viscosity of the liquid, a is the particle 
diameter and t>j is is the velocity of particle i. The iner- 
tia terms are neglected so this becomes a force balance 
equation, which is solved for the displacements of the 
particles. The conservative forces are the sum of a steep 
core potential, and the depletion force obtained via the 
following two body potential energy term. 



U 



(3) 



where the first term, Ujjci is the power law repulsive 
hard core, 



Uhc 
kT 



(4) 



n = 36 is used to avoid anomalies that may result from 
having a softer potential pZfl while retaining computer 
time efficiency. The two-body depletion potential used 



was: 



= QH(L, nj 



-A 



(•5) 



Where r*y is the separation between spheres i and j. If 
k is the size (diameter) ratio between a polymer and a 
particle at the polymer volume fraction <f> p , then 



II. METHODS 

A. The model 

The Brownian Dynamics model is based on a Langevin 
type equation of motion for each colloid particle with 
a fluctuating random force accounting for the thermal 
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Strictly speaking, the parameter <fi p = ^irk 2 o~ 2 nR, where 
ur is the polymer number density in a reservoir of pure 
polymer which would match the chemical potential of the 
polymer solution with colloid particles at volume fraction 
<j) H. It useful that the potential used here is of finite 
range, as this avoids truncation errors and allows us to 
define nearest neighbours and interacting particles pre- 
cisely. 

The Brownian forces represent the action of the liquid 
molecules averaged out over one time step. They average 
out to zero, and are uncorrelated in time. 



it can be measured in light scattering experiments [|20| 
It is calculated from simulation data by 



f.(?,M») = ^I] e, ' , ' (n(t ' ) " r,(t ' +t)1 ' ( 10 ) 

i 

where t w is the waiting time for that sample. The wave 
vector q is discretised to correctly take into account the 
periodic boundary conditions. Therefore, for each axis 
we have: 



< F B > = (7) 
<F B (t)F B (t')> = 2nak B TS(t-t') (8) 

The fluctuation-dissipation theorem relates the magni- 
tudes of F B and a by enforcing that every squared sep- 
arable degree of freedom satisfies the equipartition theo- 
rem. 

In order to simplify the calculation reduced units are 
used in the simulation, k B T becoming the energy unit. 
The polymer volume fraction, <j) p , is the only control pa- 
rameter for the depth of potential, and the ratio of the 
colloid to polymer size k controls the range of the poten- 
tial. This was fixed at k = 0.1 for all simulations shown 
here. 

All distances are measured in units of the particle di- 
ameter a. 

The natural time unit to use is the Brownian relaxation 
time. It is defined as the time for a particle to diffuse its 
own length scale by Brownian motion. 
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with L the linear dimension of the box. To improve the 
statistics the results are averaged over the x, y and z axes. 
Since it depends both on t and t w , Fg is a so-called two 
time quantity, which are most useful in determining and 
characterising any aging behaviour if it is present. 

In order to calculate particle correlations in real space 
the radial distribution function g(r) is calculated, as well 
as the self part of the Van Hove correlation function 
which samples particle displacements in time. 



1 N 

(r) = ^-^^^r^-r^-r) 



where < n > is the average number density, 



For a micron size particle in water, tr is of order 0.2 
seconds. Note however the strong dependence on parti- 
cle size. The typical length of runs showed here was of 
5000 tr and the longest run was twice that. This is a 
great increase compared to previous studies which were 
typically run for about 100 tr. 

The system being simulated is a three dimensional cu- 
bic box with 4000 particles in Periodic Boundary Condi- 
tions. 



B. Quantities of interest 

As mentioned earlier the self part of the intermediate 
scattering function (also called incoherent intermediate 
scattering function) gives information about the relax- 
ation of density fluctuations in the sample. It lets us 
assess whether particles are displaced on a length scale 
set by 2l or whether kinetic processes of that wavelength 
are 'frozen out', it is a commonly calculated quantity as 



G s (r,t) = -Y,<S(\n(t)-n(0)\-r)> 

i 

where < ... > is an average over time origins. The quan- 
tity 



S(r,t) = 4nr 2 G s {r,t) 

gives the probability of finding a particle at a distance r 
after a time t, and G s (r, 0) = g(r). 
In the Gaussian approximation, 

GsM)= (I^)i ""{wi) 

which is a solution of the diffusion equation. 

Our results are divided into two parts; the first part 
will be concerned with 'macroscopic' quantities, i.e. cor- 
relation functions and quantities that have been aver- 
aged over several systems; the second part will look more 
closely at the 'microscopies' of individual systems. 
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III. RESULTS 
A. Macroscopics 

As in previous simulations 0, [l5| the samples were 
equilibrated at high temperature (<f>= 0) and instanta- 
neously quenched to an appropriate value of <j) p . Work in 
|L5j determined that for <j) p > 0.3 crystals were no longer 
obtained. For these relatively deep quenches, amorphous, 
non-equilibrium structures were obtained; these networks 
come under the general name of gels or transient gels 
when seen experimentally before they collapse gravita- 
tionally. They are qualitatively similar to created ex- 
perimentally gels and LJ gels from simulations. All gels 
considered here have <f> = 0.30 and (j> p = 0.40, giving a 
potential minimum of U m i n = —SkT. Properties are av- 
eraged over five samples unless otherwise stated. 

In figure |l| the potential energy of the system, U, is 
scaled by U — > n Ij U . giving the approximate number 
of neighbours assuming inter-particle distances are near 
the potential minimum. The gel particles have approxi- 
mately seven neighbours. It is immediately obvious that 
we are dealing with a non-equilibrium system since, af- 
ter the initial very fast drop in energy after quenching, 
a long regime of slow decay is entered; no equilibrium 
value was ever reached in the simulations we ran. Such 
small decays can be fitted with a number of functions, 
however for comparison with Barrat and Kob |4lf| , it was 
fitted with a power law U = a - t~ b , as shown in the inset 
of figure [j], yielding b = 0.3. 

It is revealing to look at the total number of bonds in 
the system as gelation occurs. When plotted on log-linear 
scale (figure |2[) three different growth regimes are evi- 
dent, the first up to IOt^ , another from IOt/j to around 
100ta , and the last one from there on. Note that since 
the ordinate range is small, a log-log plot shows practi- 
cally the same result, hence one cannot easily distinguish 
between a low power law growth and logarithmic growth. 
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time / x R 

FIG. 3: Beta relaxation at short times for qd — 
23.0, 32.9, 46.0, with stretched exponential fits given in text 
(white line). Inset: Mean Squared Displacement (black cir- 
cles) for the same time scale with fit to equation [ill (white 
line) . 



The fits are of the form N = c.log(t) +d where the values 
of c are shown on the plot. As explained in the discus- 
sion we are here mostly concerned with behaviour over 
the last of these regimes. 

The ISF was found to display the expected two step 
relaxation at short ((3 relaxation) and long (a relaxation) 
times. They are treated separately here since the aver- 
aging cannot be done in the same way. For the waiting 
time independent (3 relaxation an average over time ori- 
gins can be used, whereas for the a relaxation the non- 
ergodic nature of the system means only an ensemble 
average can be performed. Figure |3| shows this initial de- 
cay for increasing values of q. As the gel forms, the mean 
squared displacements of the particles is Fickean, but as 
aggregation proceeds the behaviour becomes subdiffusive 
(< Ar(r) 2 > oc t e , e < 1). Finally, for measurements 
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20 40 60 

qd 

FIG. 4: Non-ergodicity parameter as a function of q, with 
exponential fit. 

in the region of interest the mean squared displacement 
takes on the form of figure ^ (inset) where a plateau is 
reached at displacements that are small compared to the 
particle size. This is interpreted in terms of the Krall- 
Weitz model for the initial decay of the correlation func- 
tion | p3[ . This model considers the relaxation due to 
the elastic modes present in the gel, with the particles 
assumed to be confined to the gel strands. Their con- 
clusion was that the mean squared displacement should 
obey: 

< Ar(r) 2 >=5 2 [l-e { ~^ r ] (11) 

The particle displacements are constrained to an am- 
plitude 5; the stretched exponential form comes about 
from the sum of the relaxations on all length scales. 

Consistently with the model we are using, the (3 decay 
was fitted with the following law which is closely related 
to a stretched exponential, 

F s (q,t) = exp(g 2 < Ar 2 (t) > /6) (12) 

for all q vectors, with < Ar 2 (t) > from equation [H]. 
The results are shown in figure || for three different values 
of q. The value of the exponent was determined from 
the fits to data in figure || with equations to be 

p = 0.4 ± 0.1. Similarly, S = 0.075 ± 0.004d and 773 = 
0.26 ± 0.0At r . 

S 2 being the maximum mean squared displacement 
(excluding any a processes), therefore sets the height of 
the plateau in F s via F s (q,t) = ± £. e -? 2 (Ar?)/6 = y* 
where /* is the non-ergodicity parameter for the incoher- 
ent scattering. Figure shows the non-ergodicity param- 
eter decays exponentially for increasing q; S 2 can then be 
extracted since f q = e~( qS ^ / 6 in the gaussian approxi- 
mation. The previously extracted value of 5 was used for 
the fit. However there is a small ambiguity about where 
exactly f q should be measured, since a processes are also 



1 10 100 1000 10000 

time / x H 

FIG. 5: q vector dependence of the intermediate scattering 
function for long times and fits to equation [l^. From top curve 
qd=2.30, 3.28, 4.27, 5.26, 6.25. The waiting time is 1000 tr. 
Inset collapse of the same data when scaled as described in 
text and plotted against stretched time. 




1 10 100 1000 10000 



Time / t r 

FIG. 6: Intermediate scattering function at fixed qd = 6.25, 
and increasing waiting times, from left t w = 20, 100, 500, 
1500, 2250 tr. 



present there is no true plateau; in this case the values 
of f q were measured after 1.5t/{. 

We now turn to the long time dynamics; the second de- 
cay (alpha relaxation) is the part of the dynamics which 
is expected to show signs of aging. Because the system is 
non-ergodic, an ensemble average cannot be replaced by 
a time average. Unfortunately, very long runs are needed 
in order to get to the time scales that are relevant for this 
study, and therefore only a limited number of runs can 
be performed due to computer time limitations. The av- 
erages were therefore performed over a set of five runs. 
Information about the dynamics on different wavelengths 
can be extracted by looking at data from a fixed waiting 
time and sampling at different q. For a given q the wait- 
ing times were then varied to observe the aging in the 
dynamics. 
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Figure |5j shows the alpha relaxations for a single wait- 
ing time and the q vectors qd = 2.3,3.28,4.27,5.26, and 
6.25. The relaxations were fitted to a stretched exponen- 
tial law: 



F s (q,t) = Aexp(~((t-t w )/T a r) (13) 

where A where was a constant chosen to be the plateau 
height and /i was kept constant for different q. The data 
collapses onto a master curve when plotted as a function 
of stretched time (t — t w Y once normalised by A and r q 
(inset of figure ^|). The relaxation times T a i p h a (q) which 
are extracted from the fits for different q where found to 
follow a power law relationship 

T a ~ q~ V 

with v = v(t w ). The best fit values extracted were 
slightly age dependent as summarised in Table |[ 
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TABLE I: Summary of the exponents as a function of age of 
the gel 

For waiting times t w < 1000 tr, as well as fast de- 
creasing values of /i, it was noticed that the quality of 
the fits decreased becoming impossible for t w < 500 tr. 
The a relaxation at these early times displayed a rela- 
tively fast decay followed by a slower one (see curves for 
t w = 20, 100 tr in figure ||. These relaxations were fitted 
with a double stretched exponential: 

+A M:(^rf)} (14) 

The justification for this law is given in the discussion 
section. Once again the relaxation times contain wave- 
length dependent information, r' a2 ~ q" 2 ; from the fits 
v' 2 ~ 2.2 for all appropriate waiting times as summarised 
in table H. r' a2 was found to depend only very weakly on 
q. The exponent fi' 2 = 0-7 was held constant and f q was 
used to give the correct plateau height. 
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FIG. 7: Intermediate scattering curves made to coincide at 
F s = 0.71. 

The significance of these exponents is discussed in the 
next section. 

Figure ^ shows the aging behaviour of the gels. The 
data was all collected at a single q value, qd = 6.25. The 
smallest waiting time considered here is 20 tr. As the 
waiting times were increased from 20 tr up to 2250 tr, 
a lengthening of the (3 plateau was observed, with the 
relaxation of the system becoming slower with increasing 
waiting time. This is classic aging behaviour for non- 
equilibrium system and glasses and has been seen for 
supercooled LJ systems by Barat and Kob 

The collapse of the a decay is often observed in exper- 
iments on glasses where one often finds that each set of 
data can be scaled using a single time scale which encom- 
passes the whole parameter dependence [p2| . Hence the 
following scaling: 

F s (q,t,t w ) ~ F(t/T 7 ) 

with t 7 = Tj(t w ). Using an arbitrary value of F s = 
0.71 to define a relaxation time for each curve t 7 = 
t(F s = 0.71), the results were rescaled and reploted on 
figure |?] using waiting times > 100 tr. The lower limit 
was chosen to be at the start of the last regime seen in 
figure ||. As can be seen the data in the range considered 
appears to collapse well. This relatively good collapse is 
surprising in that the exponent of the stretched exponen- 
tial /i varies with waiting time. 

Figure || shows the dependence of r 7 on t w . The inter- 
pretation of the results in figure || is slightly ambiguous. 
It is possible to interpret this data in two ways. Firstly in 
terms of a single power law regime r 7 ~ with £ = 1.26. 
Secondly in terms of two power law regimes, the first with 
£ = 1.33 and the second with £ = 0.62, the separation 
between the regimes being around 1000 tr. All of these 
possible values for £ are not far from 1, and similar to 
aforementioned results. It was also checked that the same 
behaviour resulted if values of T a were extracted directly 
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FIG. 8: Relaxation times extracted from the rescaling of fig- 




FIG. 9: Crystalline Order Parameter. Inset: distribution of 
particle coordinations in gel; solid line: t w = 3000 tr; dot- 
dashed line: tw = 100 tr 



from the fitted stretched exponentials. The issue of one 
or two regimes will be discussed further on. 

Another scaling was proposed by Mussel and Rieger 
for the LJ supercooled fluid results |Q and was at- 
tempted for comparison. F s (q,t,t w ) ~ F{ln[(t + 
t w )//jt]/ln(t w /n)}, where is a fit parameter and plays 
the role of an affective microscopic length scale. This 
was proposed by Fisher and Huse in the context of spin 
glasses [Q. This scaling was also applied to our data, 
however it did not produce a successful collapse. 



B. Microscopies 

Over the period of time considered, it is seen that al- 
though the average numbers of bonds in the system in- 
creases, there is little increased crystalinity as shown by 
the local order parameter, denoted ip x (main part of fig- 
ure 0). The local crystaline order parameter is defined 
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FIG. 10: Self part of the Van Hove correlation function. Main 
graph: long time interval, t w = 3000 tr. Inset: short time 
interval, tw — 10 tr. The peak area is fitted with a gaussian, 
shown as a dot-dashed line. 
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FIG. 11: Displacement correlation function, with an expo- 
nential fit. 

as the fraction of of particles which are 12 coordinated 
or neighbouring a 12 coordinated site. The inset of fig- 
ure |] shows the probability density distribution of the 
coordination of the particles in the system. The initial 
crystalinity, which is present very soon after the quench, 
is just the result of the high initial colloid volume frac- 
tion together with the strong local order, which is itself 
due to the lack of angular rigidity of the bonds. There is 
therefore no ambiguity about a gel with a non-zero crys- 
talline local order parameter. At later times an increase 
of the average position of the peak is seen but no sig- 
nificant increase in crystalinity. This already rules out 
one of the possible hypothesises for the aging of the gel, 
namely slow growth of crystals around the crystallites 
present from formation. 

In order to obtain a better understanding of the re- 
laxation mechanisms at work the particle dynamics are 
sampled in real space through the Van Hove correlation 
function. S(r,t) gives the probability of finding a particle 
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at time t + St at a distance r from its position at time 
t. Due to the non-ergodicity it is not possible to time 
average this correlation function either. However the ac- 
curacy of this result was checked by averaging over a time 
interval that was short compared by the typical a relax- 
ation time scale. Figure |l^ show a pair of typical results; 
they were calculated for t w = 3000 tr and St = 50 tr 
and 3000 tr. For these parameters we expect to be in 
the (3 plateau for the former time interval and in the a 
relaxation regime for the latter. Indeed, the peaks of the 
distributions are around 0.05d and 0.16d which are con- 
sistent with earlier deductions about those regimes. Also 
shown on figure [l0| are the Gaussian fits for the peaks 
and the regions around them. The difference in the be- 
haviour of the dynamics in those two regimes can here 
be seen more intuitively. For times up to and includ- 
ing the (3 plateau, the gel particles' displacements are 
nearly Gaussian, consistent with a diffusion-like mecha- 
nism. The main part of the plot shows the displacements 
during the a relaxation regime, the distribution is very 
much broader with the longer displacements not fitting 
within a Gaussian. The displacements found after the 
peak follow a power law S(r,t) ~ r - ^, with / = 2.8 ±0.4. 

The strongly non-Gaussian displacements rule out a 
diffusive mechanism, and hint at cooperative motion. 
This can be seen visually when colour coding particle 
displacements on a 3D representation of the gel (not 
shown here). Clusters of particles are observed to be 
displaced by similar amounts, in fact this was also true 
when colour coding displacements which were less then 
the peak value. In order to show the presence of collec- 
tive displacements clearly, a velocity correlation function 
was calculated: 
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FIG. 12: Close up the intermediate scattering function; qd ■ 
6.25, t w = 500 tr. Inset: same for a 500 particles system. 
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Cv{r) = — J^<S(r- (n - r 3 ))v{n) ■ vfa) >, (15) 

ij 

Where W — < vfc) ■ v (rj ) > is used to normalise 
C v . This is shown in figure [ll] for t w = 4000r^ and ve- 
locities calculated over an interval of At = 50 tr. The 
data are fitted with a an exponential allowing a char- 
acteristic correlation length r v = l.9d to be extracted. 
An interesting point was that the results for the velocity 
correlation function, for all times beyond the (3 relax- 
ation, were practically unchanged, thereby showing that 
whatever the relaxation mechanism, the same degree of 
cooperativity is present. 

Now it is known that particle displacements are spa- 
tialy correlated, their temporal distribution can be stud- 
ied. For this we calculated the intermediate scattering 
function for a single gel, with no averaging procedure. In 
figure |l2| the (3 relaxation and the plateau are not ap- 
parent due to the linear time scale which is useful for 
comparison with other microscopic results. It is appar- 
ent that the relaxation processes at work are not evenly 
distributed in time. A succession of plateaus and sharp 
drops are observed. Although the displacements are co- 



FIG. 13: Close up of the scaled total number of bonds in the 
system and a 10 point running average. 



operative that this function is averaged over all particles 
in the system, therefore a sharp drop in the averaged 
function is a sign of major cooperative rearrangement. 
For a smaller sized system where a greater fraction of 
the particles is included in any given strand the plateaus 
are closer to the horizontal and the fast decorrelations are 
sharper (see inset of figure [l^ for a 500 particle system). 

Following on from the localisation in time and in space 
of cooperative motion, we found similar non-uniform be- 
haviour for the total number of bonds in the system, a 
close up of which is shown on figure [l3] (scaled to show 
the average number of neighbours). A running average 
has been used to better show the underlying trends of 
the data. Here again fast increases are interspaced with 
plateaus and short time decreases. Qualitatively we can 
see that both the size and gradient of the rapid increases 
tend to decrease as time goes on. This was also true of 
F s and the same behaviour is apparent in closeups of the 
energy (not shown here). It is consistent with a slowing 
down of relaxation processes as the gel ages. 
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Considering that the particles are strongly localised, 
it was surprising to see substantial decreases in the to- 
tal number of bonds in the system. The lifetime of a 
bond based on a simple diffusion over a barrier argument 

3100 tr for a pair of parti- 
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cles interacting with this potential. When taking into 
account the average number of neighbours within inter- 
action range, one would not expect to find many broken 
bonds in any given time interval. Due to the finite range 
of the potential, it is possible to define bonds exactly as 
being within the interaction distance; therefore, we used 
a procedure that kept track of all particles that have 
had broken or created bonds within a given time inter- 
val. Remarkably, this shows that a high fraction of the 
particles break and create bonds during the evolution of 
the system. Nine percent of the particles were found to 
have broken bonds (slightly higher for the created bonds) 
when comparing configuration files 500 tr apart, this is 
much higher than what would be expected from figure 
|l3| . This number is probably an underestimate since it 
does not include bonds broken and recreated between the 
same pair of particles over that time interval. Over any 
sizable time interval there is always a surplus of created 
bonds. In order to illustrate the effect of the dynam- 
ics we show on figure |l^ the proportion of particles that 
change their near neighbour environments compared to a 
reference configuration chosen every 1000 tr. There are 
three points of importance. There is a short time interval 
where about 5 — 7% of the particles change environment 
within a couple of tr. This regime is independent of 
the reference configuration and corresponds to the (3 de- 
cay. Its independence suggests that there is a continuous 
process of breaking and creating bonds at short times. 
This will be further discussed in the next section. Fig- 
ure [l4| also has a second regime which is waiting time 
dependent and mostly follows the trends of the aging 
behaviour, the rate at which new particles change envi- 
ronment slowing down with increasing waiting time. The 
third point to notice is the absolute values. For the top 
curve the changes in the gel involve 27% of the particles, 
a huge number considering the relatively small distances 
moved by the particles. 

The presence of particle displacements that are greater 
than the contribution from the elastic modes, of corre- 
lated dynamics and a high fraction of particles changing 
environments, meant we expected to find some correla- 
tions in the position of those particles which have bro- 
ken bonds or changed environment. However breaking 
bonds in a concentrated system is always going to lead to 
some correlation between the particles; to see this imag- 
ine pulling away one of the particles of a tetrahedron, the 
three particles left behind have all had a broken bond 
but are highly correlated in space. A visual inspection of 
the gel with highlighted particles where bonds had been 
broken confirmed the suspected correlation; we also cal- 
culated a pair correlation function <7f,(r) by selecting the 
particles that had broken bonds in a given time interval 
(figure 15), and comparing it to the normal pair correla- 
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FIG. 14: Fraction of particles which change their nearest 
neighbour environment for increasing waiting times. 
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FIG. 15: Radial Distribution function done with all the all 
the particles in the gel and by selecting particles with broken 
bonds. 



tion g(r). The correlation shows up in <?&(r) as a higher 
first neighbour peak and a large peak beyond it which is 
found at the limit of the interaction potential, which is 
no surprise since it is by definition the distance at which 
one considers a bond to be broken. In terms of the previ- 
ously given argument about pulling apart a tetrahedron, 
the higher first peak is due to the strong correlation of 
the particles left behind, while the second peak is due 
to the particles that are just beyond the interaction dis- 
tance, i.e. the ones that have just broken off. This second 
peak at the limit of the potential is probably related to 
the initialy rapid increase seen in figure since it most 
likely shows that a lot of particles are hovering near the 
edge of the potential well, accounting for much of the 
bond breaking as they pop in and out when the strand 
is deformed by thermal motion. 

Exactly the same observations were made about the 
particles that had created bonds in any reasonable time 
interval. 
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FIG. 16: Averaged normal stress tensor components (negative 
of pressure); raw data in black, 20 point running average in 
white. 

Since it was the case in other studies ||(|, one pos- 
sible driving force for the aging of the gel could be the 
internal stresses. It is known that gels formed in Periodic 
Boundary Condition have non-zero normal stresses p5[ , 
and experimentaly a gel which has been detached from 
the walls of its container shrinks macroscopicly j2Qj; fig- 
ure [l6] shows the average of the normal stresses, raw data 
in black and a running average in white. The raw data 
show the high frequency noise due to the brownian fluctu- 
ations, but the presence of underlying low frequency fluc- 
tuations is revealed by the running average. The general 
trend is that of a decrease, and there are large fluctua- 
tions in the average, although not as large as any of the 
fluctuations in the raw data. A spatial stress correlation 
function, C a {R) —< a(r)a(r + R) > was calculated but 
was found to show only very low spatial correlations. 

The distribution of bond lengths in the gel was found 
to be quasi-equilibrium with only a small discrepancy at 
distances beyond the potential minima due to the normal 
stresses present in the system, as can be seen on figure 

B 

IV. DISCUSSION AND CONCLUSIONS 

This study of the dynamics of gels formed from col- 
loids interacting via a depletion potential has showed that 
complex behaviour arises during the aging. 

It was found that the total potential energy of the sys- 
tem could be described in terms of a power law with an 
exponent b — 0.3 which is an intermediate value to what 
was obtained by Barrat and Kob for a Leonard Jones 
glass (b=0.14) II and by Parisi @ (b=0.7) for a MC 
generated glass. However, the total number of bonds in 
the system was found to grow as three separate logarith- 
mic regimes. Some of this three regime structure was also 
present in the energy plot demonstrating the difficulty 
in describing accurately these very slow growths/decays. 
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FIG. 17: Distribution of bond lengths in gel (circles). Black 
line is P — exp(— Uij/kT). Recall that range of the potential 
is Lid. 



The impossibility in calculating these quantities in ex- 
periments make verification difficult if not impossible. 
However, whether interpreting the slow growth of the 
number of bonds as a power or logarithmic law, the pres- 
ence of three distinct regimes is of much interest and 
demonstrates the complex behaviour associated with the 
formation and consolidation of the gel. We know from 
previous work that the time to percolate is less than 1 tr 
pTj . Following this one would expect some fast consoli- 
dation of the network during which thin strands coarsen. 
Since diffusion processes occur, by definition, on the time 
scale of 1 tr, one would expect complex and fast rear- 
ranging at relatively short times and short length scales 
as was seen by Haw et al. for 2D colloid simulations 
where compactification was obvious at short times [E7j. 
However, given that we know that the end product of 
this is not at thermodynamic equilibrium, these fast re- 
arrangement processes must be inactive at intermediate 
to long times (otherwise the equilibrium phases would 
result). Although all the regimes seen on figure |l| occur 
after percolation, it is beyond the scope of this paper to 
explain all three, our attention is here focused on the last 
one for several reasons: it is expected that the behaviour 
maybe less complex in this region; these time scales are 
more accessible to experiments thereby providing inter- 
esting comparisons on the simulation results and insights 
that may be tested on real systems. In order to study the 
aging and not the formation/consolidation behaviour the 
work showed here focuses on the behaviour in the regime 
from 100 tr after quenching onwards. 

As an aside it is worth mentioning that the noise in 
the energy of the gel goes as FT(U) 2 « l// 1 ' 4 , this value 
being closer to the 1/f flicker noise than the more usual 
l// 2 white noise. This has been linked in the literature to 
Self Organised Criticality and to cooperative behaviour 
in complex systems. However this needs to be studied 
further before conclusions can be drawn from it. 

The self part of the ISF calculated for this sys- 
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tern displayed the classic two step relaxation. Both of 
these relaxations were found to be well approximated by 
stretched exponentials with a plateau between them. 

The age invariant short time dynamics, corresponding 
to the beta relaxation was interpreted in terms of the 
Krall Weitz model for the decorrelation due to the elastic 
modes present in a gel. Although they first used it for 
describing low volume fraction fractal gels, the physical 
reality of the model is still present even though the gel 
is not a fractal object. This model has been shown to 
work well in describing gels at higher ip by Romer et al. 
p4j with reduced values of S 2 and r c reflecting on the 
compactness of the gels. 

The mean square displacements and the beta relax- 
ation were fitted using equations [ll] and 12; the value 
of the exponent was determined to be p = 0.4 ± 0.1 
which is distinctly smaller than Krall and Weitz (p = 0.7) 
for their low (f> polystyrene gel ]2^| . Romer et al. also 
found p — 0.7 for a high </> polystyrene gel p4j . How- 
ever Solomon and Varadan found p — 0.5 for a low (j) 
thermoreversable gel of silica particles [^| . In the Krall 
Weitz model p is related to the length scale dependent 
spring constant for fractal clusters via k(s) ~ and 
p = -j^i- From computer simulation f3 ~ 3.1 for DLCA 
p8| ; from Krall and Weitz, and Romer et al. p = 0.7, 
(3 = 2.3; whereas from Solomon and Varadan p = 0.5, 
therefore /3 = 1.0 for thermoreversable gels which the 
authors believed had some angular rigidity on the local 
level due to surface roughness and a different aggregation 
mechanism leading to a locally denser gel. Furthermore, 
previous studies by one of the authors on 2D gel net- 
works where trimers provided some angular rigidity |l6| , 
showed that for such systems ~ 1.1 — 1.3. All of the 
previous suggests, as proposed by Solomon and Varadan, 
the the elasticity exponent (3 is strongly dependent on the 
angular rigidity of the local network. Here p = 0.4 there- 
fore (5 = 0.66 suggesting a structure that only slowly 
becomes more compliant with increasing length scale. 
There are several factors influencing this in these sim- 
ulations. Firstly there is no imposed angular rigidity for 
the colloid interaction (which is different from experi- 
mental particles which may be rough or deform under 
the strength of the interactions). Secondly the colloid 
volume fraction used here is 0=0.3, much higher than 
most experiments; and thirdly, the interaction range is 
very small (0.1<i), this has been shown by Heyes and co- 
workers |l5| to affect the structures obtained when form- 
ing the gels. Short interaction distances resulted in the 
creation of gels where the local structure was very com- 
pact, for example it can be seen as a split second peak 
in g(r), whereas when longer ranged potentials were used 
'fluffier' structures resulted. The former is the case in 
this study, and it is obvious that a thicker strand of col- 
loids has a higher angular rigidity. Furthermore, this 
third point is strengthened by the high <f>, and probably 
negates the first point since the lack of angular rigidity 
of a single bond becomes irrelevant if there are no singly 
connected particles in the strand. This explains the high 



angular rigidity of this structure and hence the low val- 
ues of p and p. An interesting point that arises out of 
this discussion is that for a given maximum depth of po- 
tential, the rigidity of the gel may be influenced strongly 
by the range of the potential. This may explain why for 
the higher volume fraction gels of Romer et al. p = 0.7, 
the same as for the Krall and Weitz 's low volume frac- 
tion gels, although what the range of the interactions 
present in the aforementioned experimental systems was 
not clear. 

This physical picture is confirmed by direct visualisa- 
tion of the gel which shows the clusters of particles mov- 
ing relative to each other on time scales of order 1 tr. 
Furthermore this is consistent with the high number of 
particles which change environment during such time in- 
tervals. This is a direct consequence of clusters in relative 
motion within a relatively stiff structure, thereby causing 
many bonds to break and reform. 

The length scale dependent plateau height as charac- 
terised by f% was found to have a very large q-width. 
This result is very similar to Puertas et al. ]30[|, who 
simulated the approach to gels in concentrated fluids near 
the glass transition. They obtained a similar broad de- 
pendence showing again the strong confinement of the 
particles due to the short range potentials. This extends 
their results to a lower <fi than they have used. A further 
distinction is that this result is at a point on the phase 
diagram deep in the gel state whereas they were in the 
fluid leading up to the gel state. 

The a relaxation for waiting times > 500 was 
found to be well described by a stretched exponential 
F s (q,t,t w ) — A exp(—{{t — t w )/T a )^) where the exponent 
was always /i < 1.0. This type of law is often observed 
and results from the presence of relaxation modes with a 
wide distribution of time scales. This can be written as 
a Laplace transform over the distribution of relaxation 
modes Q: F(t) = p(7)e~ 7 *d7. The exponent of the 
stretched exponential is strongly dependent on the shape 
of the distribution p("f), wider distributions leading to 
lower exponents. In this study, the variation of the ex- 
ponent of the stretched exponential, /i = 0.45 — > 0.72 
as t w = 500 — > 2250t# clearly indicates a reduction of 
the width of the modes distribution. Furthermore the 
observation that the variation in fi is fast for small wait- 
ing times and levels to /i ~ 0.7 for longer waiting times 
and also the observation of the change in behaviour for 
t w < 500 where the a relaxation appeared to have a 'fast' 
decay followed by a 'slow' one, making it impossible to fit 
a stretched exponential, all point to the possibility of ^(7) 
being a bimodal distribution where the short time modes 
are widely separated from the long time modes. Hence 
the success of using a function containing two stretched 
exponentials to fit the decay of the alpha relaxation up 
to waiting times of t w = 100 tr. In this scheme each of 
the exponentials describes one part of the distribution. 
This behaviour was reproduced by numerical integration 
of the Laplace transform for distributions of modes made 
of two widely separated Gaussians. All of the correlation 
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function data can then be explained, [i and fi 2 describe 
the distribution of the 'slow' modes, /i = fi' 2 = 0.7 ± 0.2 
for all waiting times larger than the time for the slow- 
est of the 'fast' modes to relax. fi[ which describes the 
distribution of 'fast' modes increases towards 1.0 as the 
distribution of the modes that have not yet relaxed nar- 
rows. For intermediate to long waiting times, a single 
stretched exponential is recovered. 

The q dependence of the characteristic relaxation times 
T ai T aii T 'u2i extracted from the fits also contain informa- 
tion about the processes at work, v at long waiting times, 
and v' 2 = 2.2 ± 0.1; for the latter exponent this was true 
for all the waiting times back to t w = 100 tr. This is 
very close to a diffusive law, r ~ q~ 2 , yet with a con- 
sistent deviation. v[ was found to vary only very weakly 
with q highlighting the fact that different mechanisms are 
responsible for the 'fast' and 'slow' relaxations. It also 
shows that these 'fast' modes are not related to the Krall 
Weitz elastic modes which showed up in experiments as 

_2 

r ~ q . 

Comparing these relaxation processes with the much 
more studied glasses, the a relaxation behaviour is more 
complex in this case as opposed to homogeneous super- 
cooled liquids 41 . It may not be surprising that there 



should be some processes that would be specific to gels, 
considering their spatial inhomogeneity. 

This framework also explains why in experiments a 
constant value of the stretched exponential exponent is 
usually found, since the 'fast' modes observed here would 
be inaccessible. 

When considering the ISF in non-equilibrium systems, 
aging can be observed to occur as a lengthening of the 
plateau at intermediate times. The response becomes 
strongly dependent on waiting time after the quench. 
This was seen clearly in our data. By making the curves 
coincide at one point (or equivalently using the extracted 
parameters from the fits) , a collapse within a fairly short 
time window was observed. It is obvious that since the 
exponents of the stretched exponentials are slightly wait- 
ing time dependent, the collapse cannot be perfect over 
a long time scale, however it is an important feature dis- 
played by many gels and glasses Q. The dependence 
of t 7 on the the waiting time is explained from the pre- 
vious deductions about our system, taking into account 
the fact that r 7 is nearly equivalent to r a for a single 
stretched exponential describing the whole decay. The 
presence of slow and fast modes, and therefore a change 
of behaviour once the short modes are no longer present, 
was observed again in what appeared to be two different 
aging regimes in figure |^. It is in fact a confirmation of 
that change in behaviour due to the distribution of modes 
going from bimodal to single peaked. Therefore the cor- 
rect value of £ to compare to experiments is £ = 0.66. 
Previously r ~ *£, wi th £ = 0.88 glj and £ ~ 1.1 @ 
have been found for glasses, close to £ = 1 which would 
be the simple aging case. 

The experimental studies of Cipelletti et al. have re- 
ported v = 1, i.e. T a <~ g" 1 , meaning the length scale dis- 



placed scaling linearly with time. This was done for both 
attractive and repulsive systems, where it was thought in 
both cases that the relaxation of internal stresses dom- 
inate the dynamics. The theoretical argument from the 
the same studies also explained their value of the expo- 
nent \i = 1.5. The aging behaviour found was r Q ~ 
with £ = 0.9 after a period of exponential growth for the 
gels |2(J and £ = 0.78 for the repulsive system p2[ . The 
very different values of \i and v obtained here point to 
major difference the mechanisms of the aging even if the 
scaling of the relaxation times may look similar. 

The Van Hove correlation function gives us a more in- 
tuitive idea of the particle displacements in real space. 
For large observation times and long waiting times, as is 
the case in figure [l^, the system is in a regime where the 
displacements are greater than the plateau regime. For 
comparison with a system with diffusion-like displace- 
ments a Gaussian was fitted to the peak. The very large 
tail present shows just how non-Gaussian the dynamics 
are within the a relaxation regime. The displacements 
around the peak and up to 0.25c? appear to be Gaussian 
distributed, leaving the aforementioned large tail. Not all 
displacements within this tail are due to the same phys- 
ical processes. We can distinguish between two types of 
behaviour. Firstly single particle movements: particles 
that break away from the network, becoming part of the 
gas phase before rejoining the network somewhere else 
or moving along the gel by surface diffusion. These par- 
ticles make up all of the displacements greater than Id, 
but the numbers concerned are small (due to the high 
coordination of the particles in the gel), i.e. they could 
not account on their own for the a relaxation. Secondly, 
there are particles which have moved distances of around 
0.25g? < r < Id, but which never leave the network, and 
for which the displacements are correlated in time and 
space. We also note that the part of the displacements 
found after the peak follows a power law with an expo- 
nent / = 2.8 ± 0.4. This is, within errors, what was 
deduced by Cipelletti et al. and predicted by Bouchaud 
and Pitard ^9|. However, as we have just seen, their 
system is quite different being strongly attractive, and 
other parameters such as the wavevector dependence are 
different from this work. Further work is necessary before 
conclusions can be drawn from this. 

The microscopic data presented also indicates that the 
aging does not occur as continuous process but rather 
as a succession of 'events' which occur on a time scale 
which increases with increasing waiting time, as was ob- 
served on a close up of the number of bonds per particle. 
Indeed it seemed that for some short time intervals the 
number of bonds in the system decreased. These events 
involve a relatively large number of particles moving co- 
operatively since they can even be seen on a correlation 
function averaged over 4000 particles, and dominate the 
gel dynamics for time scales greater than the duration of 
the beta plateau. 

While the ideas of cooperative dynamics in non- 
equilibrium systems is not new, string-like motion of 
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chains of particles moving cooperatively having been ob- 
served in supercooled liquids |50| , the behaviour seen here 
is more reminiscent of the work of Sanyal and Sood ]5l|] 
who studied a dense binary colloid mixture in which step 
like displacements were observed. This was accompanied 
by some hop back motion as well. It is likely that a site 
hopping scenario is responsible for this, whereas back and 
forth displacements were not observed in the gel case, al- 
though such motion may be present at shorter times. 

Unsurprisingly, these correlated dynamics were linked 
to a correlation in the localisation of the particles that 
have broken or created bonds for a given time interval. 
The relatively high number of bonds broken (and cre- 
ated) in the gel when comparing configurations separated 
by a fixed time interval is an indication of the ease with 
which thermal fluctuation of the strands disrupt the net- 
work, causing a high fraction of particles to change envi- 
ronment. This may be a result of the spatially anisotropic 
nature of the gel, whereas large thermal flutuation are 
more unlikely to occur in an isotropic glassy system. 

Although the low frequency trend in the stress data 
from figure [Hj] is the result of the microscopic aging 
'events' that occur in the gel it is now thought by the 
authors to be unlikely to be the driving force. While 
it was a surprise to find that the stresses were not spa- 
tially correlated, this may point to the fact that forces 
that can be sustained by the inter-colloidal bonds are not 
very high in comparison to the stresses occasioned in the 
strands by the thermal fluctuations, even though they 
suffice to confine the individual colloid particles to the 
gel. Under the influence of thermal fluctuations strands 
can break bonds relatively easily; this is the only expla- 
nation for the the large number of bonds broken in the 
system. A further clue which supports this argument is 
the relative size of the high frequency stress fluctuations 
compared to both the average value of the stress, and the 
size of the low frequency fluctuations. In both cases the 
stresses caused by the thermal fluctuations are of much 
higher magnitude. This supports the idea that the inter- 
nal stresses cannot be responsible for the aging, and that 
the low frequency fluctuations and the general decreasing 
trend of the average stress are consequences of the aging, 
not driving factors. 

It is not yet clear why the relaxation times increase 
with waiting time, however the mechanisms invoked here 
involve breaking the network under the thermal fluctua- 
tions. Presumably, the network will break preferentially 
at the point at which it is weakest. These weak points be- 
tween the strands that make up the network would tend 
to disappear in time (they will keep breaking until they 
do). As the network become devoid of such weak points, 
we can expect the dynamics to slow down, thereby caus- 
ing the aging behaviour. Further work is underway to 
clarify this issue. 

Comparison of this work with that of Barrat and Kob, 
reveals great similarities in the overall properties of the 
relaxation dynamics in these gels and glasses, the same 
aging was seen in the correlation functions, as well as 



superposition of the a relaxation regime. However it is 
clear from this work that the microscopic mechanisms by 
which the aging occurs in the gels are different to the 
ones observed in homogeneous glasses, the dynamics in 
the latter involving cooperative motion of strings of par- 
ticles |S0) . Similarly, comparison with Puertas et al. j^o) , 
shows the same overall behaviour even though their sim- 
ulations focused on the fluid states leading up to the gel 
transition. The relaxation mechanisms (although they 
did not report any details) were very probably similar to 
what was observed in [ p0| . 

The aging behaviour predictions given by non- 
equilibrium MCT are power law regimes either side of the 
[3 plateau, a waiting dependent a relaxation time and vi- 
olation of the fluctuation dissipation theorem. Although 
the predictions are not specifically tested here, our data, 
at least qualitatively is not in disagreement with them. 
We found a stretched exponential fit to be better than a 
power law for the a relaxation but within a certain region 
near the plateau, a power law is certainly possible. 

The Potential Energy Surface approach lends itself well 
to the interpretation of our data, and comparison with 
aging glasses. In this picture, the gel would be in a regime 
in configuration space where the system samples basins 
of IS with ever decreasing energy, but with the barriers 
between Inherent Structures becoming greater than the 
thermal energy in the system. The slow hopping over 
barriers manifests itself in the gel as network rearrange- 
ments, and hence the behaviour of quantities such as the 
total number of bonds in the system. 

Based on the work of Kob et al. |^2| , who have showed 
that an effective time dependent temperature can be de- 
fined for an aging system, Sciortino and Tartaglia |3§| ] 
have proposed an out of equilibrium thermodynamics 
framework based on the IS description of the PES. In 
this picture, the free energy is separated into a basin term 
and a configuration term reflecting the separation of time 
scales in supercooled fluids. In the case of a temperature 
jump the system equilibrates to the bath temperature 
in the same basin rapidly. The exploration of configura- 
tion space takes place on a much longer time scale, hence 
the aging observed as the system lowers its configura- 
tional energy. This description is fitting for our system 
where figure [l?] shows that the system is in quasi thermal 
equilibrium while of course not being in thermodynamic 
equilibrium. For the correlation functions calculated for 
a Leonard Jones Binary Mixture, Sciortino and Tartaglia 
observed stretched exponentials for the equilibrium fluid, 
but logarithmic decays for the out of equilibrium L J fluid. 
Note that this is in fact an MCT prediction for quenches 
near an A3 point. Once again superposition of the a re- 
laxation held and relaxation time scales increasing with 
t w were found. 

The different regimes seen the total number of bonds in 
the gel (figure [j]), can then be interpreted in terms of the 



PES, much in the way of pSfland 
taneous quench the system is stil 



After the instan- 
fluid as aggregation 



begins (0 to 10 tr), this is followed by ridge dominated 
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dynamics (10 to 100 tr) and finally the aging behaviour 
of the last regime dominated by the slow sampling of IS 
of lower and lower energy. 

The similarities of this work on supercooled fluids with 
ours is striking and the only difficult point concerns the 
shape of the a relaxations. Whereas studies of the PES 
of liquids found logarithmic decays, we had two regimes 
of stretched exponentials. It is possible that is it more 
strongly influenced by the microscopic mechanisms than 
what was previously thought. Those mechanisms could 
well be different due to the localisation of the parti- 
cles within stands rather than being distributed nearly 
isotropicaly. The thermal fluctuations in strands being 
enough to break the network thereby causing the aging 
even though thermal fluctuations of a single particle arc 
very rarely enough for them to break away from a strand. 

This highlights a problem concerning general theories 
of the glass transition, if the microscopic details that 
bring about the final relaxation effect the macroscopic 
observations, as was seen when comparing our results 
with the stress driven systems of Cipelletti et al. pc[] 
and Solomon et al. |25j ] and the aforementioned results 
for supercooled fluids, then a theory like MCT which 
works with averaged quantities is unlikely to distinguish 
between them. On the other hand, there seems to be a 
greater chance to incorporate different microscopic relax- 
ation theories when using the PES approach in order to 
estimate the time to explore different basins, and under- 
stand its relation with the aging time. 



V. CONCLUSIONS 

To finally conclude this study, we can say that simu- 
lations now appear to be a reasonable way to study the 
short term aging of soft materials. It has allowed us to 
study the intermediate scattering function of an aging 
gel in order to compare it to experimental studies while 
at the same time to understand the microscopic mecha- 
nisms responsible for the aging. The conclusion is that 
for this weakly attractive, gel-forming, system the ag- 
ing is caused by discontinuous network rearrangements 
due to the thermal fluctuations. It was also possible to 
investigate and further strengthen the relationship be- 
tween gels and glasses. A possible extension of this study 
would concern the possibility of reproducing experimen- 
tal results for a more stress driven system to check the re- 
lationship between the a relaxation and the microscopic 
mechanisms. 



VI. ACKNOWLEDGEMENTS 

RA acknowledges the support of EPSRC and Colworth 
labs (Unilever Research) for funding; he would also like 
thank H. Lekkerkerker, W. Poon, A. Puertas, M. Fuchs, 
L. Starrs, L. Cipelletti, and members of the Polymers 
and Colloids group for many useful discussions. 



[1] W. Poon, A. Pirie, and P. Pusey, Faraday Discussions 

101, 65 (1995). 
[2] S. Asakura and F. Oosawa, J. Chem. Phys. 22, 1255 

(1954). 

[3] A. Vrij, Pure Appl. Chem. 48, 471 (1976). 

[4] A. Cast, W. Russell, and C. Hall, J. of Colloidal and 

Interface Science 96, 251 (1983). 
[5] A. Cast, W. Russell, and C. Hall, J. of Colloidal and 

Interface Science 109, 161 (1985). 
[6] H. Lekkerkerker, W. Poon, P. Pusey, A. Stroobants, and 

P. Warren, Europhys. Lett. 20, 559 (1992). 
[7] M. Dijkstra, J. Brader, and R. Evans, J. Phys.: Condens. 

Matter 11, 100079 (1999). 
[8] J. Felicity, M. Lodge, and M. Heyes, Faraday transac. 

(1995) . 

[9] J. F. M. Lodge and D. M. Heyes, J. Rheol. 43, 219 (1999). 
[10] J. Lodge and D. Heyes, Phys. Chem. Chem. Phys. 1, 
2119 (1999). 

[11] B. Bifsterbosch, M. Bos, E. Dickinson, J. van Opheusden, 
and P. Walstra, Faraday Discussions 101, 51 (1995). 

[12] E. Dickinson, J. of Colloidal and Interface Science 225, 
2 (2000). 

[13] M. Bos and J. van Opheusden, Phys. Rev. E 53, 5044 

(1996) . 

[14] K.G.Soga, J.R.Melrose, and R.C.Bali, J. Chem. Phys. 

108, 6026 (1998). 
[15] K.G.Soga, J.R.Melrose, and R.C.Bali, J. Chem. Phys. 

110, 2280 (1998). 



[16] A.H.L.West, J.R.Melrose, and R.C.Bali, Phys. Rev. E 

49, 4237 (1994). 
[17] M. Grant and W. Russel, Phys. Rev. E 47, 2606 (1993). 
[18] H. Lekkerkerker, J. Dhont, H. Verduin, C. Smits, and 

J. van Duijneveldt, Physica A 213, 18 (1995). 
[19] H. Verduin and J. Dhont, J. of Colloidal and Interface 

Science 172, 425 (1995). 
[20] L. Cipelletti, S. Manly, R.C.Ball, and D. Weitz, Phys. 

Rev. Lett. 84, 2275 (2000). 
[21] B. Abou, D. Bonn, and J. Meunier, Phys. Rev. E 64, 

021510 (2001). 

[22] L. Ramos and L. Cipelletti, Phys. Rev. Lett. 87, 245503 
(2001). 

[23] A. Krall and D. Weitz, Phys. Rev. Lett. 80, 778 (1998). 
[24] S. Romer, F. Scheffold, and P. Schurtenberger, Phys. 

Rev. Lett. 85, 4980 (2000). 
[25] M.J.Solomon and P.Varadan, Phys. Rev. E 63, 051402 

(2001). 

[26] W.Kob, Soft and fragile matter (Institute of physics, 
London, 2000). 

[27] J. Bergenholtz and M. Fuchs, Phys. Rev. E 59, 5706 
(1999). 

[28] J. Bergenholtz and M. Fuchs, J.Phys.:Condens. Matter 

11, 10171 (1999). 
[29] J. Bergenholtz, M. Fuchs, and T. Voigtmann, 

J.Phys.:Condens. Matter 12, 6575 (2000). 
[30] A.M.Puertas, M.Fuchs, and M.E.Cates, Phys. Rev. Lett. 

88, 098301 (2002). 



16 



[31] P.N.Pusey and W. Megen, Nature 320, 340 (1986). 
[32] P. Segre, V. Prasad, A. Schofield, and D. Weitz, Phys. 

Rev. Lett. 86, 6042 (2001). 
[33] J. -P. Bouchaud, L. Cugliandolo, J. Kurchan, and 

M. Mezard, Physica A 226, 243 (1996). 
[34] L. F. Cugliandolo and J. Kurchan, Phys. Rev. Lett. 71, 

173 (1993). 

[35] M. Golstein, J. Chem. Phys. 51, 3728 (1968). 

[36] J. P. Bouchaud, J. Phys. I France 2, 1705 (1992). 

[37] F. H. Stillinger, Science 267, 1935 (1995). 

[38] F. Sciortino and P. Tartaglia, J. Phys.: Condens. Matter 

13, 9127 (2001). 
[39] W. Poon, L. Starrs, S. Meeker, A. Moussai'd, R. Evans, 

P. Pusey, and M. Robins, Faraday Discuss. 112, 143 

(1999). 

[40] J. Melrose, Europhys. Lett. 19, 51 (1992). 
[41] W.Kob and J.-L.Barrat, Phys. Rev. Lett. 78, 4581 
(1997). 

[42] W.Gotze, Liquids, Freezing and Glass Transition (North- 



Holland, 1991). 

[43] H. Mussel and H.Rieger, Phys. Rev. Lett. 81, 930 (1998). 
[44] D.S.Fisher and D.A.Huse, Phys.Rev.B 38, 386 (1988). 
[45] J.F.M.Lodge and M.Heyes, J. Chem. Phys. 109, 7567 
(1998). 

[46] G. Parisi, J.Phys. A: Math. Gen. 30, L765 (1998). 

[47] M. Haw, M.Sievwright, W. Poon, and P. Pusey, Advances 

in Colloid and Interface Science 62, 1 (1995). 
[48] P.Meakin, I.Majid, S.Havlin, and H.E.Stanley, 

J.Phys.A:Math.Gen 17, L975 (1984). 
[49] J. P. Bouchaud and E. Pitard, European Physics Journal 

E 6, 231 (2001). 
[50] C.Donati, J.F.Douglas, W.Kob, S.J.Plimpton, 

P.H.Poole, and S.C.Glotzer, Phys. Rev. Lett. 80, 

2338 (1998). 

[51] S.Sanya and A.K.Sood, Phys. Rev. E 57, 908 (1998). 
[52] W. Kob, F. Sciortino, and P. Tartaglia, Europhys. Lett. 
49, 590 (2000). 



